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ABSTRACT 


Many  current  satellites  employ  on-off  thrusters  to  accomplish  attitude  control  tasks 
which  may  include  initial  acquisition,  rotational  maneuvers,  and  on-orbit  stabilization.  This 
work  shows  that  the  use  of  pulse-width  pulse-frequency  (PWPF)-modulated  thrusters 
provides  several  important  advantages  over  conventional  bang-bang  thruster  control 
methods,  including  less  thruster  activity  and  closer-to-linear  actuation.  The  PWPF 
modulator  is  implemented  in  simulations  using  the  Matrixx/Systembuild  software  package. 
Simulations  assuming  a  rigid  spacecraft  are  first  performed  to  compare  the  performance  of 
the  PWPF-modulated  thrust  controller  with  that  of  conventional  bang-bang  and  time- 
optimal  bang-bang  controllers.  The  discussion  is  then  extended  to  the  case  of  a  spacecraft 
with  structural  flexibility,  as  is  encountered  quite  often  in  three-axis  stabilized  vehicles  with 
large  fold-out  solar  arrays.  Simulations  for  comparison  of  the  controllers  are  performed 
using  the  flexible  spacecraft  dynamics  model.  The  control  loop  design  in  the  presence  of 
flexibility  and  possible  interaction  with  the  PWPF  modulator  nonlinearity  are  addressed. 
Using  a  describing  function  model  of  the  modulator,  stability  margin  with  respect  to  the 
structural  mode  limit  cycle  is  predicted.  Simulations  are  then  conducted  to  verify  the 
predicted  stability  margin. 
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I.  INTRODUCTION 


A .  REASON  FOR  RESEARCH 

This  work  was  performed  in  support  of  ongoing  research  in  flexible  spacecraft 
dynamics  and  control  being  conducted  under  the  direction  of  Dr.  Brij  Agrawal. 
Justification  is  presented  for  implementing  pulse-width  pulse-frequency  (PWPF) 
modulated  thruster  control  as  an  alternative  to  conventional  bang-bang  control  methods. 
The  goal  was  to  provide  smoother  control  for  improved  pointing  accuracy  with  less  thruster 
activation  and  reduced  excitation  of  the  elastic  modes  of  vibration.  Simulation  and  analysis 
are  presented  to  illustrate  the  improvement  afforded  by  PWPF  modulation  and  to  address 
the  stability  issues  which  arise  in  the  presence  of  nonlinear  actuation.  The  research 
provides  foundation  for  future  thruster  control  experiments  using  the  Spacecraft  Dynamics 
Laboratory  facilities. 

B.  BACKGROUND 

1 .  Attitude  Control  System  Overview 

A  satellite  attitude  control  system  (ACS)  is  employed  to  provide  vehicle 
acquisition,  rotation  maneuvers,  and  on-orbit  stabilization.  In  acquisition,  the  control 
system  acts  to  null  initial  angular  velocity  errors  and  bring  the  spacecraft  to  its  operational 
configuration.  This  situation  is  often  encountered  immediately  after  spacecraft  separation 
from  the  final  launch  vehicle  stage,  whereby  the  separation  event  has  introduced  initial 
angular  rates.  Rotational  (or  slewing)  maneuvers  are  sometimes  performed  to  reorient  the 
spacecraft  to  meet  mission  objectives,  or  to  prepare  for  translational  propulsive  maneuvers. 
On-orbit  stabilization  refers  to  the  nominal  mode  of  operation  whereby  the  control  system 
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typically  acts  in  the  regulator  mode  to  maintain  a  particular  pointing  attitude  with  respect  to 
the  earth. 

An  ACS  must  provide  stability  under  the  influence  of  various  types  of  external 
disturbance  torques.  These  may  include  the  secular  disturbance  due  to  thrust  misalignment 
during  a  translational  propulsive  maneuver,  or  cyclical  environmental  torques,  e.g.  from 
solar  pressure  or  gravity  gradient.  Internal  disturbances  are  applied  by  the  motion  of 
satellite  components  such  as  steerable  antennas,  sensor  suites  or  sun-tracking  solar  arrays. 

2 .  Introduction  to  Thruster  Attitude  Control 

Most  modem  satellite  attitude  control  systems  incorporate  several  types  of 
actuators  to  include  both  internal  momentum  exchange  devices  (e.g.  momentum  wheels, 
reaction  wheels,  control  moment  gyros)  along  with  a  complement  of  thrusters.  Thrusters 
are  used  in  situations  when  disturbance  torques  exceed  the  control  authority  of  the 
momentum  exchange  devices,  and  are  often  capable  of  much  faster  reorientation 
maneuvers.  PWPF-modulated  thrusters  have  been  used  to  maintain  stability  during 
translational  maneuvers  [Reference  1]. 

Momentum  exchange  type  actuators  can  usually  be  treated  as  proportional 
devices.  On  the  other  hand,  thrusters  are  typically  on-off  devices,  only  capable  of 
providing  a  single  fixed  thrust  level.  Because  this  is  inherently  nonlinear  actuation,  thruster 
control  systems  cannot  be  analyzed  in  the  same  systematic  and  straightforward  manner  as  a 
linear  system.  Response  characteristics  are  not  scalable  as  with  a  linear  system,  but  rather 
are  highly  dependent  on  input  amplitude.  By  considering  the  system  operation  for  the 
various  attitude  control  modes  (i.e.,  various  levels  of  error  signal  amplitude),  a  general 
understanding  of  the  control  system  behavior  can  be  obtained. 

The  main  concern  in  design  of  a  thruster  control  system  is  in  performing  the 
control  operations  within  the  required  accuracy,  with  the  least  amount  of  thruster  activity. 
Lower  thruster  activity  directly  translates  into  lower  propellant  consumption  and  less 
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thruster/  valve  wear.  Thus  mission  benefits  such  as  longer  mission  life  and/or  more 
available  payload  weight  for  the  same  amount  of  fuel  can  be  realized. 

3.  The  Complication  of  Flexibility 

Structural  flexibility  must  be  addressed  in  ACS  design.  In  the  presence  of  the 
flexibility,  the  possibility  exists  that  natural  frequencies  of  the  structure  can  fall  within  the 
control  bandwidth.  This  can  lead  to  control-structure  interaction  problems  which 
compromise  the  stability  and  pointing  accuracy  objectives  (and  possibly  harm  the 
structure).  The  dominant  natural  frequencies  of  vibration  must  be  predicted,  which  often 
requires  finite  element  modeling  of  the  structure.  The  nonlinearity  of  the  thruster- actuated 
control  system  can  interact  with  the  flexibility,  resulting  in  self-sustained  oscillation  of  the 
structure  and  back-and-forth  thruster  activity. 
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H.  THRUSTER  ATTITUDE  CONTROL  SYSTEM  DESCRIPTION 


A.  SIMPLE  BANG-BANG  CONTROL  SYSTEM 

The  block  diagram  representation  of  a  simple  satellite  bang-bang  control  system  with 
angular  position  and  rate  feedback  is  shown  in  Figure  1.  The  satellite,  depicted  in 
Figure  2,  is  modeled  as  a  rigid  rotational  inertia  which  is  controlled  by  a  complement  of 
thrusters  which  provide  either  a  positive  or  negative  torque  of  constant  magnitude  U  or  -  U 
[adapted  from  Reference  2], 


controller  &  thrusters  plant  (rigid  satellite  dynamics) 


Figure  1.  Bang-Bang  Attitude  Control  System 


Here  s  is  the  Laplace  variable,  /  is  the 
total  spacecraft  inertia,  6  is  the  controlled 
satellite  pointing  angle,  and  the  factor 
1/Is2  is  the  familiar  s-domain 
representation  of  the  rigid  body 
dynamics.  The  control  torque  of 
magnitude  +/-U  is  commanded  based  on 
the  sign  of  the  error  signal,  e(t),  which  is 
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a  combination  of  gained  velocity  and  position  error.  For  this  system  the  controller  is 
defined  by 


\U  sgn(e(t))  if  |c(t)|  >  a 
{ 0  if \e{t)\  <  a 


where 

e(t)  =  Bref  -  6(t)  -  adit) 


(1) 


and  a  is  the  deadband  half-width 


A  phase  plane  trajectory  for  the  control  system  of  Figure  1,  ignoring  the  deadband  for 
now,  is  shown  in  Figure  3.  Since  rate  feedback  introduces  damping  to  the  system,  the 
phase  plane  trajectory  shows  that  the  time  response  decays  toward  the  origin  where  both 
rate  and  position  are  zero.  In  the  absence  of  rate  feedback  the  system  is  marginally  stable, 
with  undamped  oscillation  occurring  with  a  unique  amplitude  for  each  unique  initial 
condition. 


Figure  3.  Phase  Trajectory  for  Damped  Thruster  ACS  with  No  Deadband 
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Note  that  the  switching  line  with  no  deadband  passes  through  the  origin  of  the  phase 
space.  Each  period  of  constant  thrust,  from  starting  point  A  to  switching  point  B,  and  from 
point  B  to  switching  point  C,  etc.,  forms  a  parabola  in  the  phase  space  [Reference  3], 

1 .  The  Limit  Cycle 

A  limit  cycle  can  be  described  as  a  self-sustained  oscillation  condition  which 
arises  in  nonlinear  systems.  Viewed  in  terms  of  the  phase  plane  trajectory,  a  limit  cycle  is 
a  closed  path  which  is  approached  from  a  starting  condition  either  from  inside  the  closed 
path  (usually  with  the  exception  of  the  origin),  or  from  the  outside.  The  system  trajectory 
will  continue  to  trace  the  limit  cycle  path  ad  infinitum  (given  that  it  is  a  stable  limit  cycle). 
In  thruster  control,  the  "rigid  body"  limit  cycle  is  due  to  the  minimum  on-time  of  the 
thruster  (the  "structural  mode"  limit  cycle  due  to  the  interaction  of  the  nonlinearity  with  a 
flexible  structure  is  discussed  in  Chapter  V) .  Consider  the  phase  trajectory  of  Figure  4, 
the  trajectory  near  the  origin  of  the  phase  space  for  the  bang-bang  controller  of  Figure  1. 
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The  deadband  affords  periods  of  coasting  which  are  represented  as  horizontal 
lines  in  the  phase  plot.  Each  vertical  segment  (actually  a  small  section  of  a  parabola) 
corresponds  to  a  thruster  pulse.  Due  to  the  minimum  on-time  of  the  thrust,  the  minimum 
impulse  and  hence  minimum  change  in  rate,  Ad^,  produces  the  ladder-step  approach 
toward  the  origin.  It  is  this  minimum  on-time  or  minimum  impulse  of  the  thrusters  which 
produces  the  limit  cycle  seen  near  the  origin.  Note  that  a  trajectory  starting  from  inside  the 
limit  cycle  boundary  will  also  eventually  end  up  in  the  limit  cycle  condition  (assuming  a 
nonzero  rate).  The  frequency  of  the  limit  cycle  is  of  great  concern  to  the  thruster  control 
system  designer  because  it  dictates  the  fuel  usage  for  the  nominal  mode  of  operation  which 
is  the  predominant  mode  through  the  life  of  the  satellite.  However,  the  desire  for  a 
minimum  frequency  limit  cycle  is  often  in  conflict  with  the  pointing  requirements  since  a 
tighter  requirement  on  6  for  instance  will  result  in  less  coast  period  and  hence  higher  limit 
cycle  frequency.  Note  that  in  Figure  3,  after  switching  point  C  the  absence  of  the 
deadband  and  associated  coast  periods  cause  "chatter"  in  the  approach  to  the  origin.  This  is 
wasteful  of  fuel,  and  upon  arrival  in  the  region  of  the  origin,  the  limit  cycle  is  of  much 
greater  frequency  than  in  the  system  with  the  deadband. 

2.  Time-Optimal  Bang-Bang  Control 

Time-optimal  bang-bang  control  for  a  rest-to-rest  maneuver  is  accomplished  using 
the  control  law 


u(t )  =  -Usgn\ 


w 


e 


o  Ocmd  +  2(7 


(2) 


where  all  terms  are  defined  as  in  Figure  1. 

For  any  initial  condition  the  phase  plane  trajectory  consists  of  two  parabolas  with  a 
single  switching  point,  corresponding  to  full-on  acceleration  until  the  half-time  of  the 
maneuver,  followed  by  full-on  deceleration  along  the  second  parabolic  trajectory  straight  to 
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the  origin.  This  maneuver  has  been  studied  and  modified  in  many  reports,  since  it  is  the 
most  time-efficient  approach. 

B.  PULSE- WIDTH  PULSE-FREQUENCY  MODULATION 

Most  thrusters  used  for  attitude  control  are  designed  for  on-off  operation.  Pulse-width 
pulse-frequency  (PWPF)  modulation  provides  a  means  of  producing  a  variable  average 
thrust  using  a  simple  on-off  thruster,  hence  an  approximation  to  proportional  control  is 
achieved. 

A  PWPF  modulator  is  depicted  in  Figure  5.  The  average  output  of  the  pulse  train  is 
approximately  proportional  to  the  demanded  torque  input.  Other  forms  of  pulse 
modulation,  for  instance  pulse  width  (or  pulse  frequency)  or  derived  rate  modulation  are 
used  in  other  applications  which  are  not  considered  here. 


The  modulator  includes  a  relay  with  dead- zone  and  hysteresis,  or  Schmidt  trigger,  and 
an  integrator  (filter).  A  limiter  after  the  filter  was  included  in  the  digital  implementation. 
The  feedback  from  the  output  of  the  relay  is  subtracted  from  the  input  signal  A  2  .  The 
formation  of  pulses  can  be  explained  with  the  aid  of  Figure  6,  which  depicts  the  output  of 
the  discretized  modulator  and  the  filter  used  in  this  study,  for  a  constant  input.  The  filter 
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integrates  the  error  signal  (A2  -A3 ).  Before  the  formation  of  a  pulse  the  feedback  from  the 
relay  is  zero.  So  the  magnitude  of  A2  determines  the  rate  of  integration,  and  thus  the  rate 
of  growth  of  the  output  of  the  filter.  When  the  filter  output  reaches  the  threshold  eon  of  the 
Schmidt  trigger,  the  relay  switches  on  and  turns  on  the  thruster.  The  feedback  signal  then 
resets  the  integrator  at  a  rate  proportional  to  ( Um-A2 )  until  the  filter  output  falls  below  e0ff 
at  which  time  the  thruster  signal  is  turned  off  and  the  feedback  signal  is  again  zeroed  out. 
Note  that  the  appearance  of  the  pulses  as  being  other  than  square  is  due  to  the  integration 
time  step  of  10  milliseconds.  The  first  pulse  is  commanded  on  at  .06  seconds  and 
commanded  off  at  .08  seconds. 


O  .OS  .1  .15  3.  J2S  J3  JS  .4  .45  .5 

ttane^ec. 


Figure  6.  PWPF  Modulator  Output  for  Constant  Input 
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Reference  1  notes  the  modulator  tends  to  contribute  control  loop  gain  reduction  at  low 

frequency  and  phase  lag  at  higher  frequency,  and  presents  design  parameters  based  on 

tradeoffs  between  phase  lag  and  internal  deadband  size: 

Km  =4.5 

r=  0.12976 

tn 

eon  ~  0-45  (3) 

^  =0-15 
^=1-0 

The  modulator  was  discretized  for  digital  simulation.  Using  the  signal  labels  from 
Figure  5,  the  filter  output  can  be  written 

(4) 

Considering  the  signal  at  two  consecutive  time  steps  n  and  (n+1),  and  replacing  the  Laplace 
variable  with  its  time  equivalent  (a  discrete-time  first  derivative),  Equation  (4)  can  be 
reformulated  as 

(A2  -  A3)K  =  A,  +  sA.J 

=  a,  oo+  Al(”  +  ~  A'-xr  (5) 

Solving  for  the  new  value  of  the  filter  output, 

Ai(n  +  \)=A2(n  +  l)T~  A^~  Kdt  +  A(n)x(l  -f)  (6) 

where  dt  is  the  sampling  interval  of  the  digital  simulation  to  represent  the  modulator,  which 
is  chosen  to  produce  a  minimum  on-time  or  minimum  pulse  width  of  dt.  The  digital 
implementation  of  the  modulator  is  included  in  Appendix  A.  Simulation  of  the  PWPF 
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modulator  was  performed  using  a  sampling  interval  of  10  milliseconds.  Figure  7 
illustrates  operation  of  the  modulator  for  two  sinusoidal  inputs  and  a  ramp. 


Figure  7.  PWPF  Modulator  Output  for  Various  Inputs 


The  modulation  factor  or  duty  cycle  is  defined  as  the  average  output  of  the 
modulator.  Pulse  trains  for  various  constant  inputs  were  analyzed  to  verify  the  expected 
linear  relationship  between  input  and  duty  cycle.  The  pulse  trains  are  included  in 
Appendix  A.  Once  the  pulse  width  and  frequency  are  determined,  the  duty  cycle  is 
simply  determined  by: 

duty  cycle  =  pulse  width  x  pulse  frequency  (7) 

The  output  characteristics  are  tabulated  in  Table  1  and  shown  graphically  in  Figure  8. 
The  gains  computed  in  the  table  are  the  average  values  of  the  duty  cycle  for  a  given  input. 
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TABLE  1.  MODULATOR  RESPONSE  TO  CONSTANT  INPUTS 


Figure  8.  Output  Pulse  Width,  Pulse  Frequency,  and  Duty  Cycle 
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The  modulation  factor  versus  input  graph  of  Figure  8  verifies  the  nearly  proportional 
operation.  Note  that  the  modulator  has  an  internal  deadband,  whereby  a  pulse  will  not  be 
triggered  for  an  input  signal  lower  than  about  0.15.  It  is  important  to  remember  that  even 
though  the  output  relationship  is  fairly  linear,  the  relay  constitutes  a  hard  nonlinearity  and 
hence  nonlinear  analysis  must  be  employed  to  study  the  details  of  control  system 
performance. 


C.  INTRODUCTION  TO  SIMULATIONS 

Digital  simulations  were  used  for  verifying  basic  control  system  operation,  comparing 
the  performance  of  different  types  of  thruster  control,  and  exploring  system  response  with 
structural  flexibility  and  variation  of  various  parameters.  In  most  cases  simplified  block 
diagrams  are  included  in  the  text  portion  of  the  thesis  for  clarity  and  ease  of  interpretation, 
with  corresponding  Matrixx/Systembuild  block  diagrams  included  in  Appendix  B.  All 
simulations  were  performed  for  the  single  axis  control  case,  neglecting  cross-axis 
momentum  coupling  effects.  Hence  the  rigid  body  equation  of  motion  is  simply: 

Id  =  Tc  (8) 


from  which  the  Laplace  transform  yields  the  s-domain  transfer  function  representation: 


Q(s )  1 

Tc(s)  Is 2 


(9) 


where  6  is  the  controlled  pointing  angle,  /  is  the  total  moment  of  inertia,  and  Tc  the 
control  torque.  For  the  case  of  thruster  control  Tc  is  approximated  as  a  square  pulse. 

Unless  otherwise  noted,  all  block  diagram  elements  were  discretized  with  a  sampling 
interval  of  10  milliseconds.  This  enforced  the  minimum  on-time  of  the  thruster,  while 
providing  a  sample  rate  fast  enough  to  sufficiently  capture  the  dynamic  response. 
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All  simulations  were  run  using  the  parameters  of  the  flexible  spacecraft  simulator 
(FSS)  experimental  system  (but  assumed  a  rigid  body  for  the  simulations  of  Chapter  III). 
The  values  used  for  total  moment  of  inertia  and  thruster  control  torque  were  1 1.4  kg-m2 
and  0.35  N-m,  respectively.  The  FSS  was  designed  as  a  scale  version  of  the  pitch  axis  of 
a  typical  geosynchronous  momentum-biased  satellite,  so  that  results  representative  of  a 
realistic  satellite  control  system  could  be  expected.  Indeed  these  values  scale  nicely  to  the 
pitch  axis  inertia  and  control  torque  values  (440  kg-m2  and  10  N-m)  for  the  operational 
satellite  of  Reference  1. 

A  pointing  accuracy  goal  of  +/-  .1  degrees  was  used  throughout  the  study  as  a  basis 
for  comparison.  This  is  a  realistic  (if  not  lenient)  requirement  for  a  communications 
satellite  at  geosynchronous  altitude. 
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HI.  COMPARISON  OF  THRUSTER  CONTROL  METHODS  FOR  RIGID 

SPACECRAFT 


Simulations  were  run  to  demonstrate  the  performance  of  the  pulse  modulated  control 
system  and  to  compare  it  with  the  performance  of  bang-bang  control  systems.  A  block 
diagram  control  system  representation  for  rigid  body  simulations  is  shown  in  Figure  9. 
Implementation  of  the  pulse  modulator  involved  replacing  the  elements  from  Figure  5  with 
a  single  Usercode  Block  in  Systembuild  which  calls  the  FORTRAN  routine  for  the  discrete 
modulator  in  each  integration  step.  For  the  time-optimal  bang-bang  control  system  the 


CONTROLLERS 


Figure  9.  Attitude  Control  Systems  for  Rigid  Body  Simulations 
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feedback  gain  and  forward  loop  gain  are  not  applicable.  A  linear  (non-thruster)  controller 
with  position  and  rate  feedback  would  bypass  the  controllers  and  the  thruster  block.  The 
PWPF  modulator  incorporates  an  additional  gain  K  which  will  be  considered  separately 
from  the  gain  kj.  The  disturbance  torque  Tj  is  included  in  some  simulations.  The  error 
signal,  e(t ) ,  includes  the  error  between  the  commanded  angular  position  and  the  actual 
angular  position,  and  the  error  between  the  commanded  and  actual  angular  rates  where  the 
commanded  rate  is  always  implicitly  zero. 

The  linear  controller  was  used  to  size  gains  for  a  particular  desired  performance.  The 
same  gains  were  then  used  with  the  nonlinear  controllers  and  for  the  simulations  with 
flexible  dynamics,  in  order  to  form  a  basis  for  comparison.  The  closed  loop  transfer 
function  of  the  system  for  the  case  of  linear  actuation  is: 


kJJ 


2  k,T  h 

S2+~ys  +  y 


(10) 


The  equation  is  of  the  same  form  as  the  general  equation  for  a  viscously  damped 
single-degree-of-freedom  system, 


2  rTr  2  X  ' 

s  +2Ccons  +  co„ 

where  £  and  con  are  the  damping  ratio  and  natural  frequency,  respectively.  Equating 
coefficients  of  (10)  and  (1 1)  it  follows  that 

and  C  =  (12) 

For  the  linear  design  a  damping  factor  of  .84  was  chosen  for  minimal  overshoot.  The 
value  of  (Qn  can  be  thought  of  as  the  control  bandwidth,  and  was  selected  as  .06  Hz  in 
order  to  minimize  interaction  with  the  first  pole  and  zero  of  the  flexible  dynamics  in  later 
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simulations  (.12Hz  and  .15Hz,  respectively).  Based  on  these  values  Equations  (12)  yield 
the  gain  values 

ki  =  1.4  N-m/rad 
kii=  6.7  N-m-s/rad 

The  response  to  a  ten  degree  step  input  with  linear  actuation  is  shown  below  in 
Figure  10.  Graphs  are  shown  for  angular  position  in  degrees,  angular  rate  in  degrees  per 
second,  the  gained  error  signal  kje(t),  and  the  phase  plane  trajectory.  A  2%  settling  time  of 
13.6  seconds  is  achieved,  with  less  than  one  percent  overshoot. 


Figure  10.  Step  Response  for  Linear  Actuation  with  Position  and  Rate 

Feedback 
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A.  THRUSTER  CONTROL  STEP  RESPONSE  (SLEW  MANEUVER) 

For  large  values  of  the  error  signal  the  system  response  with  thrust-modulation  is 
similar  to  the  simple  bang-bang  system  because  the  PWPF  modulator  tends  to  operate 
above  its  linear  range.  This  is  demonstrated  by  comparing  Figures  11  and  12  below, 
where  a  10°  slew  maneuver  has  been  commanded  from  rest  for  bang-bang  and  modulated 
thruster  control.  Both  controllers  complete  the  maneuver  in  about  24  seconds,  slower  than 
the  linear  case  but  with  essentially  no  overshoot  as  was  desired. 


Figure  11.  Slew  Maneuver  with  Bang-Bang  Controller 


The  control  torque  pulses  change  in  frequency  for  both  controllers.  This  can  be 
understood  by  referring  back  to  Figure  4,  which  depicted  the  "ladder-step"  phase  plane 
approach  to  the  final  commanded  state.  The  coast  period  (horizontal  line  in  phase  plane) 
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between  each  pulse  increases  as  the  trajectory  approaches  0  =  0,  hence  the  pulse  frequency 
decreases.  The  trajectory  rides  the  boundary  of  the  deadband  from  just  after  one  second 
until  the  completion  of  the  maneuver. 


Figure  12.  Slew  Maneuver  with  Thrust-Modulated  Controller,  K=39 


Note  that  the  modulated  maneuver  is  slightly  smoother  in  the  phase  plane.  The 
response  is  closer  to  the  characteristically  smooth  phase  trajectory  of  the  linear  system 
(Figure  10)  because  the  PWPF  modulator  is  closer  to  approximating  linear  actuation.  The 
smoothness  of  the  trajectory  is  an  important  factor  in  the  presence  of  flexibility.  An  abrupt 
change  in  velocity  curve  from  positive  slope  to  negative  slope  indicates  an  abrupt  change  in 
acceleration  and  hence  jerk  as  well.  The  faster  the  acceleration  is  reversed,  the  more  jerk  is 
applied,  which  for  a  flexible  structure  results  in  more  stored  flexural  energy.  Whereas  the 
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bang-bang  controller  only  operates  according  to  the  sign  of  the  error  signal,  the  PWPF 
modulator  responds  to  both  the  magnitude  and  the  sign  of  the  error.  The  difference  in 
phase  trajectory  smoothness  is  shown  more  dramatically  in  simulations  where  the  feedback 
gains  and  kj  are  all  set  to  unity  (Figures  13-15).  The  resulting  damping  ratio  of  .15  (for 
the  linear  controller)  produces  excessive  overshoot  but  also  a  faster  rise  time.  Both  of  these 
effects  are  seen  to  some  extent  with  the  nonlinear  controllers,  however  the  modulated 
system  provides  a  smoother  transition  between  the  reversal  of  acceleration. 


Figure  13.  Underdamped  Slew  Maneuver  with  Linear  Controller 
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Figure  14.  Underdamped  Slew  Maneuver  with  Simple  Bang-Bang  Control 


Figure  15.  Underdamped  Slew  Maneuver  with  PWPF  Modulator,  K=30 
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Note  from  Figure  15  that  the  pulse-modulated  response  would  closer  match  the 
linear  simulation  if  the  modulator  were  operating  in  its  linear  region  over  the  maneuver.  As 
the  time  history  of  gained  error  signal  shows,  the  modulator  input  is  above  the  modulator 
saturation  level  of  1.05  (see  Figure  8)  until  about  2.5  seconds,  causing  continuous-on 
operation. 

For  the  linear  control  loop,  application  of  the  final  value  theorem  indicates  that  the 
steady  state  error  to  a  step  input  is  zero.  The  PWPF  modulator  on  the  other  hand  contains 
an  internal  dead  zone  which  results  in  a  steady  state  error  (or  defines  the  amplitude  of  the 
rigid  body  limit  cycle),  such  as  can  be  seen  in  Figure  15.  The  small  steady  state  error  is 
achieved  only  by  using  a  high  value  of  modulator  input  gain  K  (=39  for  the  previous  runs). 
The  need  for  a  substantial  value  of  input  gain  can  be  understood  by  considering  the 
magnitude  of  the  signal  prior  to  the  modulator  input  gain.  The  initial  error  signal  is  due  to 
the  commanded  input  only,  since  the  spacecraft  is  initially  at  rest: 

e(0)x^=  10°(j|p)  x  1.0  =  .174 

which  is  the  maximum  value  of  this  signal  over  the  entire  maneuver,  already  quite  small 
considering  the  modulator  dead  zone  of  approximately  +/-  0.15.  Because  the  thruster  is 
full-on  for  the  first  2.5  seconds,  the  angular  rate  increases  quickly  to  nullify  the  error 
signal: 

e(~  2. 5  sec) xkl=  (10°x  1  -  6°xl  -  4.2°/sec  xl)  x  j^=  -.003 

which  is  much  too  low  to  trigger  the  modulator.  So  a  dramatic  gain  factor  is  needed  to 
keep  the  modulator  input  signal  above  its  internal  deadband.  It  was  found,  however,  that 
the  value  of  modulator  input  gain  could  be  substantially  lower,  say  K  =10,  while  still 
achieving  the  maneuver  (Figure  16). 
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Figure  16.  Underdamped  Slew  Maneuver  with  PWPF  Modulator,  K=10 


The  rise  time  is  not  affected  appreciably  by  the  lower  input  gain,  since  the  modulator  is 
still  initially  saturated,  producing  a  period  of  constant  thrust.  The  lower  value  of  K 
however  results  in  more  steady  state  error  because  the  deadband  is  reached  sooner.  As 
with  a  linear  system,  an  increase  in  the  loop  gain  (via  the  modulator  input  gain)  decreases 
the  steady  state  error  in  the  presence  of  a  disturbance  torque;  this  is  discussed  further  in 
Chapter  IV. 

Note  that  limit  cycle  behavior  was  not  seen  after  maneuver  completion  in  the  previous 
simulations.  This  is  because  the  final  thruster  pulse  brings  the  velocity  to  zero.  With  non¬ 
zero  initial  conditions  the  limit  cycle  is  seen  at  the  completion  of  the  maneuver  (Figure 
17). 
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Figure  17.  Thrust-Modulated  Slew  Maneuver  with  Non-Zero  Initial 

Conditions 

A  ten  degree  slew  maneuver  performed  using  the  time-optimal  bang-bang  controller  is 
shown  in  Figure  18.  Although  the  maneuver  is  faster  than  the  other  types  of  thruster 
control,  the  acceleration  is  by  definition  of  the  control  law  reversed  instantaneously  at  the 
half-maneuver  time.  The  effects  of  this  with  a  flexible  structure  can  be  seen  in  the  results 
presented  in  Reference  4,  which  show  the  maximum  strain  in  the  structure  occurring  just 
after  the  half-maneuver  point.  The  stored  energy  then  produces  a  whip-like  motion  of  the 
flexible  appendage  which  contributes  to  overshoot  at  the  end  of  the  maneuver.  In  rigid 
body  simulations  without  dead  zone,  the  sampling  interval  alone  constituted  enough  of  a 
delay  to  cause  switching  instability  at  the  end  of  the  maneuver.  Hence  it  was  necessary  to 
include  dead  bands  on  position  error,  (6cmd  -  90Ut),  and  rate  60ut,  in  order  to  avoid  thruster 

chatter  at  the  sampling  rate  of  100  Hz.  Even  with  the  dead  zone.  Figure  18  shows  a  slight 
imperfection  in  the  half-time  switching.  Note  that  the  maximum  velocity  reached 
(approximately  4.5  deg/sec)  is  over  twice  the  velocity  achieved  using  the  other  controllers. 
In  many  systems  rate  limiting  is  necessary,  which  of  course  diminishes  the  intended 
advantage  of  time-optimal  switching. 
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Figure  18.  Time-Optimal  Bang-Bang  Maneuver 


Sensitivity  of  the  time-optimal  control  system  to  parameter  uncertainty  was  studied. 
Note  that  the  switching  law  (Equation  (2))  contains  the  thrust  value  and  the  total  moment  of 
inertia.  Since  these  parameters  are  generally  not  known  exactly,  the  values  used  in  the 
switching  law  are  only  a  best  estimate.  Though  a  thruster  may  be  accurately  calibrated,  the 
actual  torque  produced  is  subject  to  mounting  misalignment  tolerances  and  plume 
impingement  effects.  Also,  the  moment  of  inertia  of  the  satellite  changes  over  the  orbit 
lifetime  with  the  depletion  of  propellant.  Uncertainty  manifests  as  overshoot  (or 
undershoot)  and  subsequent  unintended  thruster  firing  (see  Figures  19  and  20).  Note 
the  rigid  body  limit  cycle  which  is  established  after  completion  of  the  maneuver. 
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Figure  19.  Time-Optima!  Bang-Bang,  20%  Overestimate  of  Moment  of  Inertia 


Figure  20.  Time-Optimal  Bang-Bang,  20%  Underestimate  of  Moment  of  Inertia 
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With  the  PWPF-modulator  uncertainty  does  not  play  as  critical  a  role  since  an  estimate 
of  the  inertia  and  the  control  torque  are  not  part  of  the  control  law  formulation. 

B.  NOMINAL  CONTROL  MODE:  RIGID  BODY  LIMIT  CYCLE 
PERFORMANCE 

The  rigid  body  limit  cycle  performance  was  determined  for  each  controller  by 
simulating  the  response  to  small  initial  conditions.  Limit  cycles  of  the  same  size  in  the 
phase  plane  correspond  to  equivalent  limit  cycle  frequencies.  With  the  goal  of  keeping 
within  +/-  0.1°  pointing  accuracy,  the  deadband  of  the  bang-bang  controller  and  the  input 
gain  for  the  modulated  controller  were  adjusted  accordingly.  In  the  case  of  the  simple  on- 
off  controller  with  deadband,  the  limit  cycle  behavior  is  shown  in  Figure  21.  From  the 
position  graph  the  limit  cycle  frequency  is  determined  to  be  about  .026  Hz. 


Figure  21.  Rigid-Body  Limit  Cycle  with  Bang-Bang  Controller 
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A  limit  cycle  of  similar  frequency  is  achieved  with  the  PWPF  modulator  (Figure  22) 
with  an  input  gain  K  of  39.  Treatment  of  the  rigid  body  limit  cycle  in  more  general 
mathematical  terms  for  any  thruster  system  can  be  found  in  Reference  5,  which  also  notes 
that  with  position  and  rate  feedback,  the  limit  cycle  is  stable. 


Figure  22.  Rigid-Body  Limit  Cycle  with  Thrust-Modulated  Controller 


The  same  limit  cycle  is  established  for  both  controllers  because  minimum-impulse 
operation  is  established,  i.e.  each  vertical  portion  of  the  limit  cycle  corresponds  to  the 
minimum  .01  second  thruster  pulse  (with  A0lrin=  .02  deg/sec).  This  is  unexpected  since 
early  literature  indicates  a  conventional  bang-bang  controller  would  not  converge  to 
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minimum  impulse  operation,  but  would  rather  form  the  higher  frequency  limit  cycle  shown 
in  Figure  23  [adapted  from  Reference  6]. 


Figure  23.  Expected  Limit  Cycles  for  Conventional  (Bang-Bang)  and 
Minimum  Impulse  (PWPF)  Systems 

Though  time-optimal  switching  would  probably  not  be  used  for  nominal  attitude 
control,  a  limit  cycle  of  comparable  size  can  be  produced  by  introducing  plant  inertia 
uncertainty  (here  -20%  was  used  as  in  Figure  20)  .  Figure  24  shows  that  minimum- 
impulse  operation  is  not  established  and  the  resulting  frequency  is  determined  to  be  about 
.18  Hz,  substantially  higher  than  the  other  controllers. 
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Figure  24.  Rigid-Body  Limit  Cycle  with  Time-Optimal  Switching 


C.  DISTURBANCE  REJECTION 

The  response  of  the  various  control  systems  to  a  constant  disturbance  torque  was 
studied.  A  constant  torque  of  .015  N-m  was  applied  in  most  cases.  This  value  was 
originally  used  in  Reference  7  to  model  the  torque  applied  due  to  cabling  wind-up  in  the 
experimental  FSS  setup.  The  value  scales  nicely  to  the  thruster  misalignment  disturbance 
torque  countered  by  the  controller  in  Reference  1.  The  assumption  of  a  constant 
disturbance  torque  is  adequate  for  the  case  of  disturbance  due  to  thruster  misalignment 
during  a  translational  maneuver.  The  response  to  stochastic  disturbances  such  as  sensor 
noise  are  not  considered. 
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A  typical  phase  portrait  for  the  thruster  systems  [Ref.  5,  Ref.  6]  converges  to  a 
disturbed  limit  cycle  of  the  shape  shown  in  Figure  25.  The  time  history  of  Figure  26 
shows  that  as  the  deadband  is  reached  due  to  the  disturbance,  a  thrust  pulse  is  commanded 
to  counter  the  disturbance.  The  plot  of  6  shows  that  the  satellite  periodically  "bounces"  off 
the  deadband.  These  plots  were  generated  using  a  more  powerful  thruster  and  are  therefore 
intended  for  illustrative  purposes  only. 


Figure  26.  Typical  Response  to  Constant  Disturbance  Torque 


Using  the  system  parameters  considered  thus  far,  the  angular  position  responses  for 
the  bang-bang  and  the  PWPF-modulated  systems  (with  K=39)  are  shown  in  Figures  27 
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and  28.  Both  results  show  thruster  activity  at  4.2  Hz,  where  the  disturbance  torque  keeps 
the  satellite  at  the  edge  of  the  deadband. 


Figure  28.  Response  of  Pulse-Modulated  Control  System  to  Constant 


Disturbance  Torque 

Note  that  the  modulated  controller  does  not  provide  the  required  pointing  accuracy  of 
.1  degrees.  With  a  linear  controller  the  steady  state  error  under  the  influence  of  a  constant 
disturbance  torque  is  inversely  proportional  to  the  controller  gain.  With  the  PWPF 
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modulator  the  gain  on  the  error  signal  can  similarly  be  increased  in  order  to  lower  the 
steady  state  error.  This  is  demonstrated  in  flexible  body  simulations  in  Chapter  IV.  The 
flexible  body  stability  problems  which  can  be  encountered  due  to  the  thruster  activity  and 
the  increased  gain  are  further  discussed  in  Chapters  IV  and  V. 

A  closer  look  at  the  gained  error  signal  and  the  resulting  pulses  (Figure  29)  for  the 
PWPF-modulated  simulation  shows  that  when  the  error  signal  reaches  the  internal 
deadband  of  the  modulator  (point  A)  a  single  pulse  is  fired  which  decreases  the  error  signal 
until  the  end  of  the  pulse  (point  B).  The  error  signal  then  grows  under  the  influence  of  the 
disturbance  until  the  deadband  is  reached  again,  triggering  another  thruster  firing  (point  C). 
The  pulse  width  of  ten  milliseconds  is  consistent  with  the  response  to  the  input  of  -.15  as 
given  in  Table  1 ,  however  the  8  Hz  pulse  rate  cannot  be  realized  since  the  input  signal  is 
not  constant,  but  rather  drops  with  each  pulse. 


Figure  29.  Error  Signal  During  Constant  Disturbance  Torque 
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IV.  FLEXIBLE  SPACECRAFT  DYNAMICS  AND  SIMULATIONS 


The  simulations  in  Chapter  HI  demonstrated  the  advantages  of  PWPF  modulation  for 
the  rigid  body  case.  The  flexible  body  dynamics  are  now  discussed  and  comparative 
simulations  presented. 

A .  FLEXIBLE  SPACECRAFT  MODELS 

It  is  often  necessary  to  include  several  vibrational  modes  of  a  flexible  structure  in  a 
spacecraft  dynamical  model  in  order  to  adequately  design  the  control  system  and  to  foresee 
possible  modes  of  interaction  between  the  control  system  and  the  structure.  It  is 
instructional  to  first  examine  the  model  of  a  spacecraft  with  one  dominant  mode  of  vibration 
about  the  controlled  axis.  Consider  the  spacecraft  with  a  flexible  solar  array  pictured  in 
Figure  30  [adapted  from  Reference  8]. 


Figure  30:  Spacecraft  With  Flexibility 
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Assuming  a  single  undamped  dominant  torsional  mode  of  vibration,  the  transfer 
function  from  control  input  to  angular  position  output  can  be  represented  as: 


G(s)  = 


6b(s)=  s2+con2 

T(s)  s2Ib(s2  +  Q„2) 


(13) 


and 


1  +  rf- 


(14) 


where  s  is  the  Laplace  transform  variable,  lb  and  Is  are  respectively  the  rigid  body  and 
solar  array  inertias  about  the  pictured  axis,  06  is  the  controlled  rotation  angle  of  the  rigid 
central  body, T  is  the  control  torque,  and  K  is  the  torsional  stiffness.  The  term  (On  in  the 
numerator  is  the  cantilever  frequency  (fixed  base  frequency)  of  the  flexible  appendage,  and 
Qn  in  the  denominator  is  the  "system  mode"  or  "free-free"  mode  frequency  of  vibration, 
whereby  the  rigid  hub  is  no  longer  fixed  but  vibrates  together  with  the  flexible  appendage. 
The  s2  term  in  the  denominator  is  attributable  to  the  rigid  body  mode  of  the  spacecraft  (zero 
frequency).  An  important  observation  to  be  made  is  that  the  complex  zeros  of  the  plant 
transfer  function  correspond  to  the  fixed-base  vibrational  frequency  of  the  flexible 
structure.  So  the  transfer  function  implies  that  while  the  flexible  appendage  may  be 
vibrating  at  its  natural  frequency,  a  control  torque  is  supplied  to  keep  the  central  body  from 
rotating  [Ref.  9].  For  the  more  general  case  of  N  modes  of  vibration,  the  transfer  function 
representation  of  the  dynamics  can  be  represented  as 


r(  ,_e(j)  i  A[s2/fi>,2+(2£/ft),)s+i] 
{S)  T(s)  /S2Af[52/Qi2+(2C,/Qi>  +  l] 


(15) 


where  I  is  the  entire  spacecraft  inertia,  ft),  is  the  ith  cantilever  mode  natural  frequency,  and 
Qi  is  the  ith  system  mode  natural  frequency. 
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The  dynamics  of  Equation  (13)  can  be  derived  easily  because  of  the  simplification  that 
a  single  vibrational  mode  and  a  known  stiffness  K  could  represent  the  solar  array 
flexibility.  A  higher  fidelity  dynamical  model  of  a  flexible  spacecraft  can  be  developed  via 
the  hybrid-coordinate  representation  of  the  system  [Ref.  10],  which  includes  both  physical 
coordinates  and  modal  coordinates.  A  finite  element  analysis  is  performed  to  determine  the 
cantilever  frequencies  and  mode  shapes  of  the  structure.  For  the  general  case  of  three 
rotational  degrees  of  freedom,  the  linearized  equations  of  motion  are  [Ref.  1],  for  the  rigid 
body: 

W  +  Hd+DTq  =  Tc  +  Td  (16) 

and  for  the  flexible  appendage: 

q  +  2^coq ’  +  +  Dd  =  0  (17) 

where  I  is  the  entire  spacecraft  inertia  matrix,  H  is  the  momentum  coupling  matrix,  D  is  the 
rigid-elastic  coupling  matrix  of  the  flexible  appendage,  6  is  the  attitude  angle  vector,  Tc 
and  Td  are  the  control  and  disturbance  torque  vectors,  q  is  the  cantilever  modal  coordinate 
vector,  and tu2  is  a  diagonal  matrix  of  squared  cantilever  modal  frequencies.  The  inclusion 
of  the  damping  term  of  Equation  (17)  is  based  on  the  assumption  of  modal  damping, 
whereby  the  damping  matrix  has  been  uncoupled  and  a  single  value  for  the  modal  damping 
ratio,  £,  is  assumed  for  all  modes. 

The  flexible  dynamics  model  used  in  this  study  was  derived  [Ref.  7,  Ref.  1 1]  for  the 
experimental  flexible  spacecraft  system  (FSS)  using  the  hybrid-coordinate  formulation. 
The  system,  pictured  in  Figure  31,  free-floats  on  air  pads  and  vibrates  in  planar  motion, 
with  translational  motion  restricted  by  the  air  bearing.  Hence  the  momentum  coupling  term 
H  (Equation  16)  does  not  apply  and  /,  0,7c  and  7^  become  scalars.  Equations  (16)  and 
(17)  then  become  [Ref.  7]: 
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(15) 


;=i 

qi+2Cimiqi  +  (Qi2qi  +  Di0  =  O 

with  the  rigid-elastic  coupling  matrix  having  reduced  to  a  vector  with  each  element  based  on 
the  modal  vector  components  in  the  rigid  body  frame  (0*  and  0f  )  and  defined  for  each 

mode  i  of  the  flexible  appendage  as: 

Di  =  Jf  (*r0,y  -  yFti  )dm  ( 1 6) 

Each  Dj  is  calculated  using  a  FORTRAN  subroutine  which  operates  on  the  cantilever  modal 
frequencies  and  mode  shapes  from  the  finite  element  analysis  output;  the  integral  of 
Equation  (16)  becomes  a  summation  over  the  number  of  discrete  subbodies  into  which  the 
flexible  body  has  been  divided,  where  .xp  and  yp  locate  each  subbody  in  the  rigid  center 
body  frame. 


THRUSTER  SYSTEM 
AIR  TANK 

MOMENTUM  WHEEL 
ASSEMBLY 

MAIN  BODY 


Figure  31.  Flexible  Spacecraft  Simulator 

In  discretizing  the  system  via  the  finite  element  method,  the  number  of  modes  was 
truncated  at  six  to  obtain  a  compromise  between  reasonable  model  accuracy  and 
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computational  feasibility.  The  model  was  placed  into  state  space  in  preparation  for  digital 
simulation  using  the  Matrixx/Systembuild  software  package.  The  state  space  representation 

of  the  system  equations  is: 


x  =  Ax  +  Bu 
y  =  Cx  +  Du 

where  the  state  vector,  x ,  is  defined  as: 

x  =  [6  qx . q6  6  qx . q6]T 


(20) 


(21) 


The  output  y  is  the  vector  of  the  current  states,  hence  C  is  a  14x14  identity  matrix  and  it  is 
assumed  that  the  feedback  values  of  angular  position  and  rate  are  measured  exactly.  The 
system  matrix  form  used  in  Matrix  x/Systembuild  is 


5  = 


"A 

B 

C 

D 

(22) 


The  system  matrix  was  discretized  for  simulations  with  dt=  .01  second.  The  complete 
mathematical  descriptions  of  the  state  matrix  A  and  the  input  matrix  B  are  given  in 
Reference  1 1  (p.  20).  The  direct  transmission  matrix  D  is  set  to  zero. 

During  initial  simulations  with  the  flexible  dynamics,  the  natural  frequencies  of 
oscillation  were  found  to  vary  somewhat  from  the  values  presented  in  Reference  7  (pages 
21  and  110).  In  order  to  find  the  actual  frequencies  present  in  the  discretized  model,  a  bode 
plot  of  the  system  was  obtained  (Figure  32).  The  DBODE  command  in  Matrixx  was 
used  and  the  normalized  frequencies  were  divided  by  the  sampling  period  to  obtain  rad/sec. 
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Figure  32.  Bode  Plot  of  Discretized  Plant  Model 


Based  on  the  bode  plot  the  bode  plot  the  cantilever  and  system  frequencies  were 
determined  and  are  given  in  Table  2. 


TABLE  2.  FSS  DISCRETE  MODEL  CANTILEVER  AND  SYSTEM 

FREQUENCIES 


Mode 

Cantilever  Frequency 
(Hz)  (rad/sec) 

System  Frequency 
(Hz)  (rad/sec) 

1 

0.121 

0.76 

0.153 

0.96 

2 

0.330 

2.07 

0.410 

2.58 

3 

2.660 

16.71 

2.730 

17.15 

4 

3.550 

22.31 

3.630 

22.81 

5 

6.140 

38.58 

6.210 

39.02 

6 

16.610 

104.36 

16.900 

106.19 
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Reference  7  determined  the  first  cantilever  natural  frequency  to  be  .139  Hz,  and  the 
first  system  frequency  .175  Hz.  The  discrepancy  in  the  frequencies  could  partly  be  due  to 
the  method  of  discretization  used.  Consider  the  following  continuous  representation  of  a 
second  order  notch  filter  [Ref.  1]: 


(23) 


Figure  33  shows  the  continuous  frequency  response  plotted  with  the  discrete-time 
response  using  two  different  discretization  methods  (with  w—  .96).  Note  the  shift  in  zero 
frequency  introduced  by  the  zero  order  hold  method,  whereas  the  first  order  hold  method 
provides  much  better  results.  The  same  shift  may  have  occurred  when  the  FSS  dynamics 
model  was  discretized.  (This  illustration  underscores  the  difficulty  in  designing  a  notch 
filter  for  lightly  damped  poles  since  only  a  small  inaccuracy  in  the  knowledge  of  the  true 
pole  location  would  result  in  poor  filter  performance.) 


(a)  <*>) 

Figure  33.  Continuos  to  Discrete-Time  Conversion  of  Notch  Filter  using 
(a)  zero  order  hold  method,  (b)  first  order  hold  method. 


The  state  space  representation  of  the  system  contains  the  information  to  form  transfer 
functions  between  the  torque  input  and  each  of  the  fourteen  outputs.  The  transfer  function 
for  the  central  body  rotation  angle  output  was  given  in  Equation  (15),  and  can  be  modified 
into  the  form 


rrc,  0(*>-  J_rrf Vt  [I+k^s+q)?] 


Using  the  values  of  the  cantilever  and  system  natural  frequencies  from  Table  2,  and 
an  assumed  modal  damping  ratio  of  .4  percent,  the  transfer  function  for  the  first  three 
modes  becomes: 


23  *  ( 5 2  +.  0065+.  58  )(s2  +.  0 1  Is  +  4. 4  l)(s2 +.  1 34s  +  279) 
(r  +.  008.V+.  922)(.v2 +.  02  Is  +  6. 76)(.$'2  + 1 38.v  +  296) 


(25) 


The  poles  and  zeros  of  the  complete  six-mode  transfer  function  were  computed  using  a 
Matlab  routine  which  is  included  in  Appendix  C  and  are  shown  in  Figure  34. 


Real  Axis 


Figure  34:  Poles  and  Zeros  of  Transfer  Function  from  T  to  6 
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Note  that  the  underdamped  modes  appear  as  alternating  poles  and  zeros  near  the 
imaginary  axis.  Reference  9  shows  that  with  position  and  rate  feedback  or  lead 
compensation,  stability  of  the  modes  is  achieved,  that  is  the  root  locus  from  each  open  loop 
pole  to  its  adjacent  zero  will  lie  in  the  left  half  plane.  With  angular  position  and  rate 
feedback  the  loop  transfer  function  is  multiplied  by  the  factor  k]( zs+1) .  The  resulting  root 
locus  plot  near  the  origin  is  shown  in  Figure  35,  for  r  =  4.77.  An  underlying  assumption 
is  that  the  sensor  and  actuator  are  colocated  at  the  "rigid"  central  hub  of  the  spacecraft. 
With  flexibility  between  the  sensor  and  actuator,  lead  compensation  can  actually  cause 
structural  mode  instability  (see  Reference  12,  Appendix  A.4). 


Figure  35.  Root  Locus  vs.  Forward  Loop  Gain  With  Position 

and  Rate  Feedback 
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If  the  controller  gain  were  chosen  to  locate  the  closed  loop  poles  near  points  A  for 
maximum  damping  of  the  flexible  modes,  the  resulting  rigid  body  control  would  be 
considerably  underdamped.  More  realistically,  the  closed  loop  poles  would  probably  be 
placed  with  the  rigid  body  pole  closer  to  point  B,  for  nearly-critical  damping  of  the  center 
body.  Hence  the  flexible  poles  would  be  dictated  by  the  rigid  body  controller,  and  would 
be  much  nearer  the  imaginary  axis.  Although  the  stable  root  locus  plot  is  guaranteed  with 
the  linear  controller,  nonlinearity  can  cause  closed  loop  pole  locations  to  change  as  a 
function  of  time.  Nonlinear  analysis  is  discussed  in  Chapter  V. 
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B .  COMPARISON  OF  THRUSTER  CONTROL  METHODS  FOR 

FLEXIBLE  SPACECRAFT 

Comparisons  of  bang-bang  and  pulse-modulated  thruster  control  in  the  presence  of 
structural  flexibility  are  presented  in  this  section.  The  first  two  of  the  six  cantilever  modal 
coordinates  are  plotted  to  give  some  feel  for  the  response  of  the  flexible  appendage  with 
respect  to  the  central  body.  The  simulations  were  run  at  the  same  sampling  rate  (100  Hz) 
and  with  the  same  basic  block  diagrams  as  the  previous  rigid  body  simulations,  only  the 
flexible  body  dynamics  were  substituted  for  the  rigid. 

1 .  Step  Response 

The  ten-degree  step  response  is  shown  in  Figures  36-41  for  the  PWPF- 
modulator ,  bang-bang,  and  time-optimal  bang-bang  controllers.  The  time  responses  of  the 
first  two  modal  coordinates  are  included  to  provide  some  indication  of  the  arm  excitation. 
The  PWPF-modulator  input  gain  of  30  (Figure  36)  was  reduced  to  10  for  the  run  of 
Figure  37.  Note  the  reduced  thruster  firing  with  little  loss  in  accuracy,  as  well  as  the 
rapid  decay  in  the  modal  coordinate  amplitudes.  The  oscillation  in  pointing  angle  is  due  to 
the  input  gain  to  the  modulator,  which  increases  the  loop  gain.  Hence  the  original  linear 
design  with  low  controller  bandwidth  to  avoid  interaction  with  the  first  mode  is  not 
realized.  It  is  possible  that  further  "optimization"  of  the  modulator  could  help  improve  this 
aspect  of  the  results. 

For  the  first  bang-bang  run  (Figure  38),  the  dead  bands  were  left  as  were 
required  to  provide  the  minimum-impulse  limit  cycle  from  the  rigid  body  simulations  (with 
+/-.10  excursion  in  6).  The  resulting  thruster  actuation  is  unsatisfactory,  although  tight 

control  of  the  center  body  is  noted.  The  forward  deadband  was  opened  up  to  +/-  .2°  to 
produce  the  result  of  Figure  39.  The  increased  deadband  reduces  the  "density"  of 
thruster  firings,  however  the  oscillations  and  firings  continue  until  almost  200  seconds. 
Similar  results  were  obtained  using  the  time- optimal  controller.  The  first  run  (Figure  40) 
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shows  a  much  faster  maneuver  (as  is  expected),  however  excessive  switching  occurs  at  the 
termination.  Again  opening  up  the  deadband  (Figure  41)  provides  virtually  no  benefit  in 
terms  of  reducing  thruster  firings,  and  results  in  undesirable  oscillations  in  pointing  angle. 
The  phase  plane  plot  is  shown  in  this  run  to  point  out  the  limit-cycle-like  behavior, 
however  further  investigation  showed  thruster  firing  ends  after  about  200  seconds.  The 
delayed  overshoot  occurring  just  after  six  seconds  is  probably  due  to  the  "whipping" 
motion  of  the  flexible  appendage  as  it  releases  the  energy  stored  during  the  maneuver.  The 
amplitude  of  cantilever  modal  coordinate  qi  reaches  .26,  much  higher  than  for  the  thrust- 
modulated  maneuver  where  the  amplitude  is  .085.  The  rate  of  decay  is  also  more  gradual 
than  the  modulator  case. 


tkna<aae) 


Figure  36.  PWPF-Modulated  Step  Response,  K=30,  Flexible  Dynamics 
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modulator  output 


Figure  37.  PWPF-Modulated  Step  Response,  K=10,  Flexible  Dynamics 


Figure  38.  Simple  Bang-Bang  Step  Response,  Flexible  Dynamics 
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Figure  39.  Simple  Bang-Bang  Response  with  Increased  Deadband 


Figure  40.  Time-Optimal  Bang-Bang  Step  Response,  Flexible  Dynamics 
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Figure  41.  Time-Optimal  Bang-Bang  Response  with  Increased  Deadband 


Though  these  simulations  were  intended  to  point  out  the  advantages  of  PWPF 
modulated  thruster  control,  it  is  important  to  remember  that  variations  of  each  control  law 
have  been  devised  to  improve  the  performance  in  the  presence  of  flexibility  [e.g.  Agrawal 
&  Bang,  Dodds].  The  basic  idea  in  slew  maneuvers  is  to  minimize  the  post-maneuver 
energy  in  the  flexible  structure,  however  this  desire  tends  to  conflict  with  the  fact  that  a  fast 
slew  is  also  desired.  [Hailey]  demonstrated  a  slewing  method  which  developed  a  more 
gradual  angular  acceleration  using  torque-shaping  and  resulted  in  less  excitation  of  the 
structure.  This  method  and  others  would  show  similar  improvement  for  any  of  the 
controllers. 
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2.  Small  Initial  Condition  Response 

Figures  42-45  show  the  responses  of  each  of  the  control  systems  to  small 
initial  conditions,  where  phase  plane  trajectories  are  included  rather  than  the  modal 
responses.  Thruster  activity  in  the  first  PWPF  case  (Figure  42)  was  reduced  further  by 
implementing  an  external  deadband  of  +/-  .1°  before  the  modulator  (Figure  43).  The 
tradeoff  is  seen  to  be  slightly  less  accuracy,  but  still  well  within  the  goal  of  +/-  .1°.  The 
simple  bang-bang  response  (Figure  44)  is  comparable  to  the  modulated  response.  The 
time-optimal  bang-bang,  which  implements  the  dead  bands  directly  in  the  feedback  loop,  is 
seen  to  bounce  back  and  forth  within  the  +/-  .1°  deadband. 


Figure  42.  PWPF-Modulated  Response  to  Small  Initial  Conditions, 

Flexible  Dynamics 
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modulator  output 


Figure  44.  Simple  Bang-Bang  Initial  Condition  Response 
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Figure  45.  Time-Optimal  Bang-Bang  Initial  Condition  Response 
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3 .  Disturbance  Rejection  with  PWPF-Modulated  Controller 

The  behavior  of  the  modulated  control  system  during  a  constant  disturbance  is 
illustrated  in  Figures  46-50.  The  disturbance  of  .015  N-m  was  applied  for  20  seconds 
(except  for  Figure  46),  starting  at  one  second.  Figure  46  shows  that  with  K- 30,  the 
required  pointing  accuracy  appears  to  be  achieved.  After  the  disturbance,  a  rigid  body  limit 
cycle  of  frequency  .05  Hz  appears  to  have  been  established.  The  longer  run  of  Figure  47 
(70  second  disturbance)  clearly  shows  that  the  required  pointing  accuracy  is  maintained. 
The  results  of  Figures  48-50  show  that  as  the  input  gain  is  increased,  the  steady-state 
error  is  reduced,  however  some  degree  of  modulator  instability  is  shown  for  the  case  of  K— 
200  (Figure  50).  After  the  disturbance,  the  over-excited  modulator  still  manages  to  drive 
the  system  toward  zero  pointing  error. 


Figure  46.  Modulated  Controller  Response  to  Disturbance  Torque,  K=30 
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ire  47.  Modulated  Controller  Response  to  Disturbance  Torque,  K=30 


Figure  48.  Modulated  Controller  Response  to  Disturbance  Torque,  K=50 
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V.  NONLINEAR  CONTROL  LOOP  STABILITY  ANALYSIS  IN  THE 
PRESENCE  OF  STRUCTURAL  FLEXIBILITY 


The  stability  of  the  control  system  is  a  function  of  all  elements  of  the  control  loop.  The 
various  elements  other  than  constant  gain  add  phase  lead  or  lag  (along  with  gain 
contribution  except  in  the  case  of  pure  time  delay),  and  hence  contribute  to  the  overall  gain 
margin  and  phase  margin.  As  was  previously  stated,  for  a  control  loop  with  linear 
actuation  and  position  and  rate  feedback,  the  stability  of  each  mode  is  assured.  Though  the 
control  loop  may  induce  oscillations  in  the  structure,  the  oscillations  will  experience 
exponential  decay.  The  discontinuous  nature  of  the  PWPF  modulator,  however,  can  cause 
the  control  to  interact  with  the  dynamics  and  result  in  self-sustained,  or  limit  cycle, 
oscillations. 

A.  LIMIT  CYCLE  DETERMINATION  USING  DESCRIBING  FUNCTION 
ANALYSIS 

In  order  to  predict  the  existence  and  characteristics  of  limit  cycles  in  a  system, 
describing  function  analysis  can  be  employed.  Consider  the  block  diagram  depiction  of  a 
system  which  has  been  separated  into  linear  and  nonlinear  elements  (Figure  51,  adapted 
from  Reference  2). 


Nonlinear  Element 

Describing  Function  Linear  Element 


r(t) 


^  x(t) 

N(A,  co) 

w(t) 

L(jco) 

y(t) 

J 

t 

Figure  51.  Control  System  with  Nonlinear  Element 
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For  a  sinusoidal  input  to  the  nonlinear  element,  i.e.  x(t)=  Asin(cot),  the  output  w(t)  of  the 
nonlinear  element  is  generally  periodic  and  hence  can  be  represented  by  a  Fourier  series 
expansion.  The  describing  function  method  makes  the  approximation  that  only  the 
components  corresponding  to  the  fundamental  harmonic  of  the  expansion,  nco  with  n=  1, 
are  to  be  considered.  This  approximation  is  based  on  the  assumption  that  the  linear 
element  has  low-pass  properties  and  hence  tends  to  attenuate  the  higher  harmonics  of  w(t). 
The  describing  function  of  the  nonlinear'  element  is  defined  as  the  complex  fundamental- 
harmonic  gain  of  a  nonlinearity  in  the  presence  of  a  driving  sinusoid  [Ref.  14]: 


N(A,co)  = 


phasor  representation  of  output  component  at  frequency  co 
phasor  representation  of  input  component  at  frequency  co 

A1  (A,  CO)  a,  (a,o>) 

A 


=  jib  +  j<0 


(26) 


where  the  magnitude  and  phase  of  the  describing  function  are 
\N(A,co)\  =  j^Ja12+b12 

A  (27) 

ZN(A,co)  =  (j>  =  tan'^Oj/fy) 

The  describing  function  N(A,<a),  is  analogous  to  the  frequency  response  function, 
H(jco),  of  a  linear  system.  But  whereas  the  linear  frequency  response  function  is 
independent  of  the  amplitude  of  the  input  signal,  the  nonlinear  "transfer  function"  or 
describing  function  depends  on  both  the  input  amplitude  and  frequency.  In  linear  control 
loop  design,  closed-loop  poles  can  be  placed  to  assure  a  stable  response.  With  a  nonlinear 
element,  however,  the  closed-loop  pole  locations  vary  with  time  due  to  the  gain  variation 
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with  input  amplitude.  Hence,  the  root  locations  in  the  complex  plane  can  oscillate  between 
locations  in  the  left  and  right-half  plane,  giving  rise  to  limit  cycle  behavior. 

By  representing  the  nonlinear  element  by  a  describing  function  approximation,  some 
of  the  same  graphical  techniques  from  linear  control  analysis  can  be  employed  to  predict  the 
limit  cycle  behavior  and  to  determine  stability  margins  with  respect  to  the  limit  cycle 
condition.  Consider  the  configuration  from  Figure  51,  assuming  that  a  self-sustained 
oscillation  of  amplitude  A  and  frequency  co  is  established  and  hence  the  input,  r( t)  =  0. 
The  following  relations  can  be  obtained  from  the  figure: 

x  =  -y 

w  =  N(A,co)x 
y  =  L(j(0)w 


and  hence 

y  =  L(jco)N(A,co)(-y) ,  or 
y(l  +  L(jco)N(A,  eo))  =  0 

Since  y*0,  the  limit  cycle  frequency  and  amplitude  must  satisfy  the  relationship 

w>"mz>  (28) 

A  bode  plot  representation  of  a  system  with  nonlinearity  is  depicted  in  Figure  52. 
The  negative  inverse  of  the  describing  function  gain  and  phase  are  superimposed  with  the 
linear  plot.  The  gain  and  phase  margins  are  redefined  with  respect  to  the  nonlinear 
boundary,  as  opposed  to  the  -180  degree  phase  and  0  dB  gain  lines.  The  limit  cycle 
frequency  and  amplitude  corresponds  to  the  point  which  satisfies  Equation  (28),  the  point 
of  intersection  between  the  nonlinear  and  linear  gain  plots. 
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Figure  52.  Gain  Margin  (Gm)  and  Phase  Margin  (Pm)  for  (a)  Linear 
System,  and  (b)  System  with  Nonlinear  Element 


B  .  ANALYSIS  OF  PWPF  MODULATOR 

For  the  rigid  body  case,  the  limit  cycle  frequency  was  easily  determined  by 
simulation.  With  a  flexible  spacecraft  the  dynamic  characteristics  of  the  modulator  need  to 
be  considered  in  order  to  predict  the  limit  cycle  behavior  which  can  result  from  the 
interaction  of  the  nonlinearity  of  the  PWPF  modulator  with  the  control  system.  The 
dynamic  characteristics  can  be  analyzed  using  the  describing  function  method. 

Although  the  describing  function  method  is  based  on  the  assumption  that  only  the 
fundamental  harmonic  of  the  nonlinear  response  is  significant,  the  residual  higher 
harmonics  have  been  considered  [Ref.  14].  The  corrected-conventional  describing  function 
uses  the  root-mean-square  value  of  the  first  and  third  harmonics,  with  the  phase  of  the 
fundamental.  A  corrected-conventional  single  input  describing  function  (SIDF)  for  the 
PWPF  modulator  is  presented  in  Reference  1,  for  which  the  gain  and  phase  of  the  SIDF 
boundary,  N(Amjco),  are  determined  to  be 
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Gain(dB)  =  20  log(5  /  Am ) 
Phase(deg )  =  -57.3tan"'(7^£u) 

where 

h-Uon-  Uoff  =  hysteresis  width 


,  h  1V(7»2-h 

on  l  +  exp(A/(7»)J  Km 


B  = 

3  = 

S3  = 

h  =  Uon~  Uoff  =  hysteresis  width 
co  -  input  frequency  (rad  /  s) 

A  =  minimum  pulse  width  (s) 


(29) 


A  block  diagram  representation  of  the  continuous-time  control  system  is  shown  in  Figure 
53. 


PWPF  Modulator  Flexible 


Figure  53.  PWPF-Modulated  Thruster  Control  Loop 
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Bode  plots  for  the  control  loop  were  developed  in  Matlab  using  the  continuous-time 
dynamics  description  from  Equation  (24).  The  code  is  included  in  Appendix  C.  The 
bode  plot  for  the  first  five  pole-zero  pairs  of  the  plant  dynamics  is  shown  in  Figure  54. 
The  gain  plot  agrees  closely  with  that  found  in  Matrixx  using  the  DBODE  command  on  the 
discrete-time  state  space  model  of  the  plant  (see  Figure  32).  The  phase  plot  from  Matrixx, 
however,  shows  phase  lag  at  the  higher  frequencies  which  is  not  seen  in  Figure  54.  This 
is  thought  to  be  another  consequence  of  the  discretization  process.  The  possible  effect  that 
this  discrepancy  could  have  is  that  the  predicted  phase  margin  with  respect  to  a  limit  cycle 
may  be  smaller  than  the  analysis  using  the  Matlab  model  would  indicate. 


Figure  54.  Bode  Plot  of  Flexible  Dynamics 
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Figure  55  shows  the  plant  with  position  and  rate  feedback  control  (with  K=30) 
superimposed  with  the  SIDF  boundary  developed  using  Equations  (29).  The  nonlinear 
boundary  intersects  the  linear  plot  for  the  first  five  modes,  with  linear  control  loop  gain 
reduction  for  the  sixth  system  mode  frequency  sufficient  enough  that  limit  cycle  interaction 
is  not  predicted.  Hence  five  possible  limit  cycle  interactions  are  indicated;  only  the  third 
structural  mode  limit  cycle  at  «  21  rad/sec  was  found  to  be  significant,  as  will  be  discussed 
shortly.  It  is  interesting  to  note  that  in  Reference  1  only  the  first  mode  pole  (representing 
the  first  system  mode  of  vibration)  needed  to  be  considered  because  the  other  modes  of 
vibration  were  far  removed  in  frequency  from  the  first  and  thus  far  below  the  SIDF 
boundary.  For  the  flexible  spacecraft  model  in  this  study,  however,  the  modal  density  was 
high  enough  that  higher  frequency  flexible  mode  limit  cycles  could  be  expected.  Hence  a 
fairly  accurate  plant  model  is  desired  since  exclusion  of  the  higher  frequency  modes  in 
initial  analysis  would  fail  to  predict  the  higher  frequency  limit  cycles. 


Figure  55.  Control  Loop  Bode  Plot  with  Nonlinear  Boundary,  K=30. 
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Note  that  the  position  and  rate  feedback  effectively  provide  lead  compensation  which  is 
indicated  in  Figure  55  by  the  increase  in  gain  of  the  higher  frequency  underdamped  pole- 
zero  pairs  (in  addition  to  the  32  dB  contributed  by  the  total  gain  factor  of  42),  and  by  the 
phase  lead.  The  phase  lead  is  desirable  as  it  provides  greater  phase  margin  with  respect  to 
the  nonlinear  boundary.  The  gain  increase  however  causes  the  linear  boundary  to  intersect 
with  the  SIDF  boundary  at  the  higher  frequencies  where  limit  cycle  interaction  is  of  even 
more  concern  since  higher  frequency  limit  cycles  are  more  wasteful  of  fuel.  With  the  lead 
effect  modeled  as  zs+1,  a  decrease  in  the  rate  feedback  gain  t  would  increase  the  break 
frequency  and  hence  lessen  the  gain  increase  seen  for  the  third-through-fifth  pole-zero 
pairs.  But  it  follows  from  Equation  12  (Chapter  III)  that  this  would  also  decrease  the 
damping  ratio,  which  may  produce  undesirable  system  response  (i.e.  too  much  overshoot). 

The  stability  of  the  control  loop  with  respect  to  the  third  structural  mode  limit  cycle  can 
be  evaluated  by  adding  pure  time  delay  to  the  control  loop,  which  introduces  phase  lag  with 
no  gain  contribution.  With  a  large  enough  time  delay  the  phase  margin  is  effectively 
eliminated  and  limit  cycling  is  predicted  to  occur.  Pure  time  delay  in  the  s-domain  is  by 
definition: 

G(s)  =  e~Ts  (30) 

where  T  is  the  delay  time. 

Employing  a  Taylor  series  expansion  where  higher  order  terms  are  neglected,  an 
approximation  to  pure  time  delay  with  sample  and  zero-order-hold  [Ref.  1]  is 

1  -  e~Ts  1 

■  .  g-s-r - - -  (31) 

Ts  7V/12  +  7s/2  +  l 

Note  that  for  this  analysis  the  use  of  continuous-time  descriptions  of  the  control  loop 
elements  is  founded  because  the  sampling  rate  used  in  the  discrete-time  simulations  is 
sufficiently  higher  than  the  dynamic  response  frequencies. 
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Figure  56.  Loop  Bode  Plot  with  60  msec  Delay,  K=30. 
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C .  SIMULATION  OF  LIMIT  CYCLE  BEHAVIOR 

The  analytical  stability  margin  was  verified  by  adding  loop  time  delay  to  the  simulation 
block  diagram  and  performing  simulations  for  various  delay  times.  In  real  satellite  control 
system  implementation,  the  on-board  processing  capability  often  limits  the  sampling  rate 
and  the  speed  with  which  each  subsequent  control  signal  command  is  computed.  A  multi¬ 
rate  sampling  scheme  was  incorporated  to  simulate  the  microprocessor  sampling  interval 
(which  contributes  to  the  delay).  The  block  diagram  of  Figure  57  incorporates  the  two 
sampling  rates  and  the  process  delay.  The  microprocessor  sampling  interval  was  chosen  to 
be  40  msec.  The  total  time  delay  is  a  combination  of  the  microprocessor  sampling  interval 
and  pure  computational  delay.  The  PWPF  modulator  continues  to  cycle  at  100  Hz 
(providing  pulse  command  updates  every  10  msec),  but  receives  an  input  signal  every  40 
msec.  Hence  the  modulator  sees  a  constant  input  for  four  10  msec  cycles  at  a  time. 


PWPF  Modulator 


Process 


Td 


Figure  57.  Higher  Fidelity  Control  Loop  Model 


The  structural  mode  limit  cycle  was  excited  in  the  simulations  first  by  placing  initial 
conditions  on  the  controlled  angle  (0)  and  the  modal  coordinate  (qi)  corresponding  to  the 
deformed  shape  of  the  satellite  and  flexible  appendage,  for  whichever  mode  was  closest  to 
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the  predicted  limit  cycle  frequency.  It  was  found,  however,  that  just  about  any  initial 
condition  would  trigger  the  limit  cycle  in  the  presence  of  the  proper  amount  of  time  delay. 
Figures  58  and  59  shows  the  time  history  of  the  physical  and  modal  coordinates, 
respectively  with  initial  conditions  on  6,  qi  and  q2  ,  and  a  60  msec  time  delay.  A  limit 
cycle  of  frequency  17.3  rad/sec  is  established  after  about  four  seconds,  close  to  the 
predicted  frequency  of  about  20  rad/sec  from  the  analysis  of  Figure  56.  The  limit  cycle 
amplitude  (from  the  q  plot)  is  of  acceptable  magnitude,  within  the  pointing  accuracy 
"requirement"  of  +/-  .1  degrees.  Note  the  decay  of  all  modal  coordinates  except  q3,  which 
grows  in  amplitude  to  the  limit  cycle  condition.  The  PWPF  modulator  output  shows 
sustained  alternating  positive  and  negative  thruster  firings,  characteristic  of  the  limit  cycle 
condition. 

The  modal  coordinate  time  histories  could  be  converted  into  physical  arm  deflections. 
However,  the  modal  magnitudes  are  somewhat  smaller  than  in  the  results  from  Reference 
7,  indicating  that  the  physical  deflections  are  not  of  concern  from  the  standpoint  of  the 
possibility  of  structural  damage.  Also,  the  main  concern  here  is  the  control  of  the  central 
body. 
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PWPF  output  error  signal  r«te{dpe) 


Figure  58.  Initial  Condition  Response  with  K=30,  60  msec  Delay 


Figure  59.  Initial  Cond.  Modal  Response  with  K=30,  60  msec  Delay 


66 


Figures  60  and  61  show  that  with  a  reduced  time  delay  of  20  milliseconds  the  limit 
cycle  is  not  excited.  Note  the  .12  Hz  oscillation  in  6  ,  the  frequency  of  the  first  cantilever 
mode.  It  was  not  possible  to  excite  limit  cycles  near  the  fourth  and  fifth  system  modes,  as 
predicted  by  the  analysis. 


Figure  60.  Initial  Condition  Response,  K=30,  20  msec  Delay 
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Figure  61.  Initial  Condition  Modal  Response,  K=30,  20  msec  Delay 


The  simulation  of  Figures  62  and  63  was  performed  with  no  time  delay.  Note  that 
after  the  thruster  firing  stops,  the  body  oscillates  at  the  first  system  mode  frequency  of  .15 
Hz.  An  important  implication  of  on-off  thruster  control  can  be  noted  from  the  modal 
coordinate  response.  Note  that  each  thruster  firing  excites  the  higher  modes,  and  the 
modes  ring  down  between  pulses.  This  points  out  one  of  the  shortcomings  of  trying  to 
think  of  an  on-off  thruster  system  in  linear  response  terms.  Even  though  the  output  in  8 
appears  to  be  typical  of  a  sinusoidally  forced  linear  system,  each  thruster  pulse  input  to  the 
plant  actually  has  high  frequency  content,  as  we  know  from  the  trigonometric  series 
expansion  of  a  square  wave.  Note  that  the  sampling  rate  of  100  Hz  is  six  times  the  highest 
system  frequency,  so  that  proper  capture  of  the  qt  dynamics  is  assured. 


68 


PWPF  output  error  signet  retype)  po«Mon(<l«g) 


Figure  62.  Initial  Condition  Response  with  K=30,  No  Delay 
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Figure  63.  Initial  Condition  Modal  Response  with  K=30,  No  Delay 
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Figures  64  and  65  show  that  a  slew  maneuver  can  still  be  performed  in  the  presence 
of  the  third  system-mode  limit  cycle. 


Figure  64.  Slew  Maneuver  with  K=30,  60  msec  Delay 


Figure  65.  Slew  Maneuver  Modal  Response,  K=30,  60  msec  Delay 
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phase(deg)  magnitude(dB) 


The  analysis  was  continued  by  considering  a  lower  PWPF  input  gain,  predicting  the 
limit  cycle  frequency  from  the  bode  plot,  then  adding  delay  until  limit  cycling  occurred. 
Figure  66  is  the  linear  bode  plot  with  nonlinear  boundary  for  the  case  K=10.  For  the 
case  with  K=30  the  fourth  and  fifth  system  mode  limit  cycles  were  not  excitable,  hence  the 
SIDF  boundary  is  overly  conservative.  Assuming  the  same  conservatism  for  the  case 
K=10,  it  was  predicted  that  only  the  second  mode  limit  cycle  would  be  evident.  The 
predicted  frequency  from  the  bode  plot  is  5.5  rad/sec. 


Figure  66.  Loop  Bode  Plot  With  Nonlinear  Boundary,  K=10,  No  Delay 
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With  a  delay  of  280  milliseconds,  the  gain  margin  with  respect  to  the  second  mode 
limit  cycle  is  eliminated  (Figure  67).  The  gain  attenuation  after  10  rad/sec  is  due  to  the 


Figure  67.  Loop  Bode  Plot  With  Nonlinear  Boundary 


K=10,  280  msec  Delay 

The  predicted  stability  margin  with  respect  to  the  second  mode  limit  cycle  was  verified 
and  the  results  are  shown  in  Figures  68  and  69.  The  limit  cycle  frequency  is  4.7  rad/sec 
(.75  Hz).  Simulations  were  run  with  delays  of  up  to  270  milliseconds,  but  limit  cycle 
behavior  was  not  established  until  280  milliseconds  of  delay  was  added.  The  assumption 
that  the  SIDF  formulation  is  conservative  was  once  again  verified  as  the  third  mode  limit 
cycle  could  not  be  excited  for  the  case  K=10.  Note  that  the  limit  cycle  amplitude  violates 
the  pointing  requirement. 
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PWPF  output  error  signal  rate(dpa)  posttlon(deg) 


Figure  68.  Initial  Condition  Response  with  K=10,  280  msec  Delay 


Figure  69.  Initial  Cond.  Modal  Response,  K=10,  280  msec  Delay 


73 


This  analysis  has  shown  that  the  control-structure  interaction  predicted  using  the  SIDF 
formula  is  overly  conservative.  The  K=30  case  predicts  limit  cycling  near  the  fourth  and 
fifth  structural  modes,  and  likewise  with  K=10  a  third-mode  interaction  is  still  predicted. 
Since  the  describing  function  method  is  by  definition  an  approximate  method  of  predicting 
instability,  this  is  acceptable.  But  the  value  of  backing  up  the  analysis  with  simulation  is 
clearly  demonstrated. 
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VL  CONCLUSIONS 


A  pulse-width  pulse-frequency  modulator  has  been  successfully  implemented  and 
shown  to  provide  pseudo-linear  operation  over  a  certain  range  of  input.  Simulations  using 
the  modulator  indicate  that  this  technique  provides  closer-to-linear  actuation  than 
conventional  bang-bang  thruster  control  methods.  This  facilitates  a  more  tailorable 
response,  which  means  that  a  single  set  of  thrusters  can  be  used  to  meet  various  control 
objectives.  The  resulting  smoother  control  results  in  less  thruster  firing  for  identical 
maneuvers  as  compared  to  simple  bang-bang  and  time-optimal  bang-bang  controllers. 

Analysis  of  the  control  loop  was  shown  to  be  important  because  of  the  interaction 
between  the  modulator  nonlinearity  with  the  flexible  dynamics.  Nonlinear  describing 
function  analysis  proved  to  be  useful  in  predicting  flexible-mode  limit  cycle  behavior. 

A .  RECOMMENDATIONS  FOR  FURTHER  STUDY 

•  The  PWPF  modulator  in  this  study  required  a  fairly  large  input  gain  to  provide 
adequate  pointing  accuracy.  The  resulting  loop  gain  caused  the  bandwidth  of  the  controller 
to  overlap  with  flexible  modes.  The  various  elements  of  the  modulator  could  be  studied 
and  further  optimized  for  the  specific  system  in  hand  in  order  to  reduce  the  excitation  to  the 
structure  introduced  during  maneuvers.  Filtering  of  the  first  flexible  mode  (as  in  Reference 
1)  could  also  be  applied. 

•  Operation  of  the  PWPF  modulator  and  comparisons  to  the  performance  of  other 
controllers  should  be  verified  by  experiments  with  the  Flexible  Spacecraft  System. 

•  Use  of  thruster  pulse  modulation  for  momentum  dumping  should  be  investigated  with 
the  goal  of  reducing  the  angular  offset  noted  in  Reference  7. 
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•  A  combined-controller  approach  should  be  investigated.  For  instance,  the  minimum 
time  solution  could  be  implemented  to  provide  initial  slewing,  with  switching  to  a  PWPF- 
modulated  controller  to  avoid  excessive  thruster  activity  at  the  termination  of  the  maneuver. 

•  The  combination  of  center  body  control  with  active  control  of  the  flexible  appendage 
using  piezoelectric  actuators  and  sensors  is  recommended.  The  "smart  structures" 
experiments  to  date  have  been  performed  with  the  appendage  cantilevered  from  a  fixed 
center  body.  A  multiple-input  multiple-output  controller  could  be  developed  which 
simultaneously  considers  end  point  control  of  the  flexible  appendage  with  the  center  body 
pointing  control. 
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InpuUl.OS 


APPENDIX  A 


The  following  graphs  depict  the  modulator  output  for  various  constant  inputs. 
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The  following  is  a  listing  of  the  digital  implementation  of  the  PWPF  modulator.  The 
code  was  implemented  within  the  format  of  the  general  Userblock  subroutine 
(USROl.FOR).  Refer  to  Figure  5  (Chapter  2)  for  definition  of  the  parameters.  The 
Userblock  is  implemented  in  simulations  by  first  compiling  the  progam  ($  FOR  USR01) 
followed  by  linking  ($  MWSLINK  USR01)  which  creates  an  appended  MWS  executable 
file  with  default  name  MYMWS.  Entrance  into  Systembuild  is  then  accomplished  using  the 
command  $MWS  MYMWS.  See  Chapter  5  of  the  Matrixx/Systembuild  Modeling  and 
Simulation  manual  for  more  information  on  Usercode  Blocks. 


c  INITIALIZATION 

IF(TIME.EQ.O)THEN 

c  Modulator  parameters 
K=4.5 
Tau=.  12976 
lm=.54 
um=l 
uon=  .45 
uoff=  .15 
dt=  .01 

c  limiter 

if(al.gt.lm)  al=  lm 
if(al.lt.-lm)  al=  -lm 

c  +fire  (schmidt  trigger) 
if(al.gt.uon)then 
ml7=  1. 
else 

ml  7=0. 
end  if 

if(al.gt.uoff)then 
nl7=l. 
else 
nl7=0. 
end  if 

if((ml7+nl7).gt.O.)  then 
ol7=l. 
else 
ol7=0. 
end  if 
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c  -fire  (schmidt  trigger) 
if(al.lt.-uon)then 
pl7=l. 
else 
pl7=0. 
end  if 

if(al.lt.-uoff)then 
ql7=l. 
else 
ql7=0. 
end  if 

if((pl7+ql7).gt.0.)then 
rl7=l. 
else 
rl7=0. 
end  if 

c  initial  output 

a3=  um*(ol7-rl7) 

y(2)=  a3  //  y  is  output  vector  of  Usercode  Block 
go  to  100 
END  IF 

c  MAIN  LOOP 

c  new  input  u  from  simulation 
a2=  u(l) 

c  filter 

al=  dt/Tau*K*(a2-a3)  +  alold*(l.-dt/Tau) 
c  limiter 

if(al.gt.lm)  al=lm 
if(al.lt.-Im)  al=-lm 

c  +  fire  (schmidt  trigger) 
if(al.gt.uon)  then 
ml8=l. 
else 

ml  8=0. 
end  if 

if(al.gt.uoff)then 
nl8=l. 
else 
nl8=0. 
end  if 

if((ml8+nl8*ol7).gt.0.)then 
ol8=l. 
else 
ol8=0. 
end  if 

c  -  fire  (schmidt  trigger) 
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if(al.lt.-uon)then 
pl8=l. 
else 
pl8=0. 
end  if 

if(al.lt.-uoff)then 
ql8=l. 
else 
ql8=0. 
end  if 

if((pl8+ql8*rl7).gt.0.)then 
rl8=l. 
else 
rl8=0. 
end  if 

c  OUTPUT 

a3=  um*(ol8-rl8) 
y(l)=al 
y(2)=  a3 

c  reassign  the  fire  commands  and  filter  output  for  next  time  step 
ol7=  0I8  F 

rl7=  rl8 
alold=al 

100  RETURN 
end 
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The  following  listing  of  a  sample  simulation/  plotting  routine  is  included  for  the  benefit 
of  the  next  Matrixx  user.  These  routines  help  minimize  the  amount  of  time  spent  retyping/ 
modifying  commands  on  the  Matrixx  command  line.  The  program  is  accessed  using 
o\\edt  filename,  and  executed  using  oexecute  ('filename').  In  the  editor  use  Ctrl  z  to 
get  back  to  the  editor  command  line,  'quit'  to  exit  without  saving  revisions,  'exit'  to  exit 
and  save  revisions,  'cti  to  get  back  to  editing. 

dt=.01; 

tmax=100; 

t=  [0:dt:tmax]';  //  Define  time  vector 
U=  l ;  //  lower  limit  for  graphs 

ul=  6000;  //  upper  limit  for  graphs 

u=  ones(t)*10*pi/180;  //  input  for  simulation  (10°  step) 

y=  sim(t,  u);  //  simulation  command 

plot(t(ll : ul) ,y (ll:ul,4) ,’upper  left,  ylabeyposition(deg)/  xlabel/time(sec)/') 
plot(t(ll : ul),y (11 : ul, 3) ,'upper  right,  ylabehrate(dps)/  xlabehtime(sec)/  date') 
plot(y(:,4),y(:,3),'lower  right,  xlabel/posittion(deg)/  ylabel/ang  vel/... 
title/phase  plane  trajectory/') 
pause 

hard  //hardcopy  of  plot 
Wprint/queue=ln03r  matplotps 

plot([t(ll:ul)  t01:ul)  t01:ul)],  ty01:ul,9)  y(ll:uM0)  y(l!:ul,l  1)]....  //strip  chart  plot 
'strip  ylabel/q  1  Iq2lq3/  xlabel/time(sec)/  ymin-1.2  ymaxl.2') 
pause 
hard 

Wprint/queue=ln03r  matplot.ps 
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Discrete  Super-Block  Sampling  Interval  First  Sample  Ext. Inputs  Ext. Outputs  Enable 

PWPF  modulated  Ctrl  multirate _ O.QlQQ  0 .  1  12  Parent 


APPENDIX  B 
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[simple  bang-bang  with  deadband] 
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Discrete  Super-Block  Sampling  Interval  First  Sample  Ext. Inputs  Ext. Outputs  Enable 

pwpf  early  implementation  0.0156  0.  1  4  Parent 


Attempts  were  made  (unsuccessfully)  to  implement  the  PWPF  modulator 
using  a  Usercode  Block  for  the  Schmidt  trigger  and  Systembuild  elements 
for  the  remainder  of  the  loop. 
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APPENDIX  C 


These  Matlab  programs  were  used  for  development  of  the  loop  transfer  function  and 
single  input  describing  function  (SIDF)  boundary  for  analysis  in  the  frequency  domain. 

%  The  zeros  of  the  plant  transfer  function  correspond  to  the 
%  eigenvalues  of  the  appendage  only,  the  cantilever  frequecies  (rad/sec). 
z(l)=  .76;  z(2)=  2.1;  z(3)=  16.7; 
z(4)=  22.3;  z(5)=  38.6;  z(6)=  104.4; 

%  Construct  the  plant  transfer  function  numerator  polynomial: 
zeta=  .004;  %  .4  percent  modal  damping  assumed 
numplant=  [1  2*zeta*z(l)  z(l)A2]; 
denfactor=  z(l)A2; 
for i=2:6 

numplant=  conv(numplant,[l  2*zeta*z(i)  z(i)A2]); 
denfactor=  denfactor*z(i)A2; 
end 

%  The  poles  of  the  plant  transfer  function  correspond  to  the  free-free  or 
%  system  modes  (rad/sec): 
p(l)=  .96;  p(2)=2.6;  p(3)=  17.2; 
p(4)=  22.8;  p(5)=  39;  p(6)=  106; 


%  Construct  the  plant  transfer  function  denominator  polynomial: 
denplant=  [1  2*zeta*p(l)  p(l)A2]; 
numfactor=  p(l)A2; 
for  i=2:6 

denplant=  conv(denplant,[l  0  p(i)A2]); 
numfactor=  numfactor*p(i)A2; 
end 
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%  Factor  in  the  rigid  body  dynamics  and  the  product-of-polesA2  to 
%  product-of-zerosA2  ratio: 

Iest=  1 1.4;  %  estimate  of  the  total  MOI 
denplant=  Iest*conv(  denplant,[  conv([l  0],[1  0])  ] ); 
factor=  numfactor/denfactor; 
numplant=  numplant*factor; 

%  Poles  and  zeros  of  transfer  function  (no  control) 
axis([-10  10-150  150]) 

rlocus(numplant,denplant)  %title('  Poles  and  zeros  of  plant  xfer  function') 

pause 

clg 

axis 


%  Bode  plot  of  plant 
w=  logspace(-l,2,400); 
bode(numplant,denplant,w),  pause 


%  Now  factor  in  the  controller/feedback  loop  to  form  the  loop  xfer  function 


K=  10; 

k=  K*  1.3985; 
tau=  4.7746; 
numH=  [tau  1]; 
denH=  1; 


%  PWPF  modulator  input  gain 
%  forward  loop  gain 
%  factor  for  rate  gain 
%  H  is  feedback  loop 


numloop=  k*conv(numplant,numH);  %  combine  plant  and  controller 
denloop=  conv(denplant,denH); 

bode(numloop,denloop,w)  %  Bode  of  loop  xfer  funtion,  GH 
pause 

[maglinear,phaselinear,w]=  bode(numloop,denloop,w);  %linear  mag  and  phase 
%[Gm,Pm,Wcg,Wcp]=margin(maglinear,phaselinear,w);  %linear  margins 


%  Determine  the  SIDF  boundary  (magnitude  and  phase): 
[magsidf,phasesidf]=  sidf;  %  calls  function  sidf 
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magsidf=  magsidf+1.5;  %  corrected  to  reflect  [Bayloc]  experimental  results 


%  Plot  both  linear  mag  and  phase  and  nonlinear  boundary: 
subplot(211),semilogx(w,magsidf,'--',w,  20*logl0(maglinear)),grid 
%title('Linear  magnitude  with  SIDF  boundary  NO  DELAY  24feb') 
%gtext('Nonlinear  Boundary')  %gtext('Gm') 
xlabel('frequency(rad/sec)'),ylabel('magnitude(dB)') 
subplot(212),semilogx(w,phasesidf,'-',w,phaselinear),grid 
%title('Linear  phase  with  SIDF  boundary’) 

%gtext('Nonlinear  Boundary') 
gtext('Pm’) 

xlabel('frequency(rad/sec)'),ylabel('phase(deg)') 

pause 

%  The  pure  time  delay  causes  a  phase  lag  which  shrinks  the  nonlinear  phase 
%  margin.  Use  sample  and  zero-order-hold  formula  from  [Wie  &  Plescia] 

T=  .280;  %  60  msec,  delay 
numdelay=[l]; 
dendelay=  [TA2/12  T/2  1]; 

%  augment  numloop  and  denloop  and  calculate  new  magnitude  and  phase: 

numloop=conv(numloop,numdelay); 

denloop=conv(denloop,dendelay); 

[maglinear,phaselinear,w]=bode(numloop,denloop,w);  %linear  mag  and  phase 

%  Plot  both  linear  mag  and  phase  and  nonlinear  boundary  again  for  comparison: 
clg 

subplot(21  l),semilogx(w,magsidf 20*logl0(maglinear)),grid 

gtext('Nonlinear  Boundary') 

xlabel('freq  (rad/sec)’), ylabel(’magnitude  (dB)') 

subplot(212),semiIogx(w,phasesidf,'-’,w,phaselinear),grid 

gtext('Nonlinear  Boundary') 

gtext('Pm  ') 

xlabel('freq  (rad/sec)'), ylabel('phase  (deg)') 


90 


SIDF  subroutine: 


function  [magsidf,phasesidf]=  sidf 
%  Modulator  characteristics  and  SIDF  analysis  from  [Wie] 

%  Static  characteristics: 

Uon=  .45; 

Uoff=  .15; 
h=  Uon-Uoff; 

del=.01;  %minimum  pulse  width 
Um=  1; 

Tm=.  12976; 

Km=  4.5; 

%  Calculate  SIDF  boundary  gain  and  phase: 
w=logspace(- 1 ,2,400); 

Am=  ( Uon-h./(  1  +  exp(del/Tm./w) )  ).*(  (  sqrt(Tm.*w)  ).A2  +  1  )/Km; 
Bl=  4/pi*Um*sin(del.*w/2); 

B3=  4/3/pi*Um*sin(3*del.*w/2); 

B=  sqrt(B  1 .  A2+B3.A2); 

magsidf=  -20.*logl0(B./Am); 

phasesidf=  -180-(  -57.3*atan(Tm.*w) ); 
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